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ABSTRACT 



<D 

Simple simulations suggest that the phase space structure of haloes identified in 
cosmological calculations is invariant under the dynamics induced by sinking sub- 
structure satellites — the background expands so as to leave the total distribution 
unchanged. We use a Fokker-Planck formulation to show that there are long lived 
solutions for densities p ~ r 7 and — 2 < 7< — 1; indices between —1 and —1.5 cor- 
' responding to the inner cusps of cosmological haloes, where the coupling is strongest; 

| steeper ones to intermediate radii. We recover the exact solutions found by Evans & 

Collett (1997); reinterpret them in terms of well defined background-satellite interac- 
CN tion; and show that these, and all other solutions, are valid for any mass spectrum 

of substructure, because the governing equation is linear in their mass weighed phase 
5^ ■ space distribution. If the spatial distribution of substructure has a milder cusp than 

the total, the system expands; when the background has a milder cusp there is com- 
. pression. It is not possible for the individual distributions to retain their original form: 

i' light particles are driven out of low energy states, being replaced by the sinking mas- 

1~ \ sive ones. If the clumps are considered solid, this takes the form of an exponential 

^**| ■ instability, with characteristic timescale of the order of the dynamical friction time, 

leading to a low energy cutoff in the distribution function of the background and a 
^ , constant density core. We show that there are long lived solutions with such a cutoff. 

They would correspond to a situation whereas the clumps are made of dense baryonic 
material. When stripping is important, as in the case of dissipationless substructure, 
it is likely that this situation is reversed — the cutoff is now in the clump distribution 
jJJ ■ function. A simple description suggests that this renders equilibria even more long 

lived. In all cases it is possible to find solutions that are long lived from the thermo- 
dynamical (energy transfer) perspective. In systems without stripping the only truly 
stable solutions however are isothermal spheres, but there are double power law solu- 
tions that may be relevant if stripping is involved. The results in this paper suggest 
that halo profiles similar to those found in dissipationless cosmological simulations are 
approximately invariant under the interaction induced by the presence of substructure 
satellites — a necessary condition for the observed 'universality'. In addition, the total 
profile, including baryons, should also be invariant; provided the latter are initially in 
the form of dense clumps, whose distribution follows that of the dark matter. 

Key words: dark matter - galaxies: haloes - diffusion - galaxies: general - galaxies: 
formation - galaxies: structure 



1 INTRODUCTION 

Collapsed structures identified within the context of dissi- 
pationless cosmological simulations are invariably found to 
exhibit some invariant, or universal, form for their density 
profiles — irrespective of their mass or the epoch at which 
they are identified (see, e.g., Navarro et al. 2004 and the 



references therein) — which, moreover, seems to reflect an 
invariant underlying phase space density distribution (Tay- 
lor & Navarro 2001). These results are potentially of fun- 
damental importance; they may signal the existence of a 
generic tendency in initially cold and nearly homogeneous 
gravitational configurations to evolve towards definite final 
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states, with a consistency reminiscent of that characterising 
the approach to a Gibbs distributions in laboratory systems, 
irrespective of the details of the initial conditions or the in- 
termediate dynamics. 

The existence of invariant final states in systems where 
the microscopic dynamics is known to be time reversible 
(Newtonian equations in our case) inevitably echoes the 
protracted debate pertaining to irreversibility in laboratory 
statistical mechanics, the main conundrum being that such 
'attractor' like behaviour is characteristic of dissipative sys- 
tems. Nevertheless, as is well known, dissipative dynamics 
can be effectively introduced in otherwise reversible systems 
by incorporating a degree of randomness. This can either be 
achieved by the presence of a fluctuating force, representing 
otherwise intractable interactions, or because the intrinsic 
dynamics is so complicated that tracking its detailed evo- 
lution from any given initial condition becomes impossible 
— the accompanying loss of information has consequences 
similar to the incorporation of random forces. Systems where 
such processes are in action may exhibit the required macro- 
scopic irreversibility that leads to invariant final states. 

Even if, in the cosmological context, the 'universality' 
may be only approximate, the relatively quaint state of af- 
fairs that transpires, rare in astrophysical research, has gen- 
erated intense interest. Discussions dealing with the origin 
of the identified density distributions can be broadly divided 
into two main theses along basic nature (the profiles result 
from the initial collapse: Lokas & Hoffman 2001; Nusser 
2001; Hiotelis 2002) or nurture (they arise from cumula- 
tive effect repeated mergers or interactions: Syer & White 
1998; Dekel Devor & Hertzoni 2003; Dekel et al. 2003) 
lines. These situations encompass the two contexts alluded 
to above: irreversible evolution can either arise from the 'vi- 
olent' internal evolution involving complicated dynamics of 
the mean field during the initial collapse or merger events; 
or phenomena arising from random encounters between, and 
dynamical friction forces on, component clumps. It is the 
latter category that concerns us here. 

The fact that haloes are always found to have basically 
the same forms for their profiles suggests that both argu- 
ments may be relevant: for, while numerical calculations sim- 
ulating the cold collapse of isolated gravitational systems do 
indeed show that the process can account for halo profiles 
(e.g., Hiotelis 2002), a halo's tenure in the context of a cos- 
mological environment will be punctuated by repeated merg- 
ers; it will contain substructure, which will be constantly re- 
plenished by continual infall. The density distribution needs 
to be invariant under the action of such interactions if the 
profile is to survive. Evidently, this will be the case if the 
underlying phase space structure is unaffected by the inter- 
actions. 

In this paper we endevour to show that this is in fact 
true; that the presence of substructure leaves the total phase 
space mass distribution of haloes with density distributions 
similar to those of cosmological-simulation-haloes (at least 
approximately) invariant. What needs to be shown is that 
when dynamical friction acts on clumps (representing dark 
matter substructure) inside a larger halo, causing them to 
sink towards the centre, the increased density that results 



from the mass deposition is precisely balanced by an outflow 
in the parent halo material; and that mass stripping from 
clumps does not change that basic conclusion. 

Since the evolution necessarily involves macroscopically 
irreversible behaviour, due to dynamical friction and random 
clump-clump encounters, the formulation presented in this 
paper may also shed light on the process of attainment of 
the universal profile via substructure interaction. Though 
this issue is only touched upon briefly here. 

As in laboratory statistical mechanics, effective ran- 
domness, born of the presence of granularity in a system, 
can be introduced by considering a macroscopic transport 
equation that contains a collision term where this random- 
ness is incorporated. In the context of gravitational systems 
the usual approximation involves the Fokker-Planck equa- 
tion, extensively used in the study of globular clusters for 
example (Spitzer 1987) . It has also been invoked in the con- 
text of cosmological haloes (Evans & Collett 1997; Wein- 
berg 2001; Ma & Bertschinger 2003). The basic assumption 
in deriving it from more basic kinetic equations is that the 
encounters are weak. This implies that it cannot describe 
the dynamics leading to major mergers between a few large 
components of roughly equal mass. The assumption however 
is well founded when what is involved is the gravitational 
interaction of a large number of clumps (representing sub- 
structure) and their motion though a smooth background 
medium (representing a parent halo) . This is the picture we 
will have in mind in this paper. 

In Section |2] we discuss simple numerical experiments 
illustrating the phase space dynamics of such a clump- 
background system. Following this, in Section |3] we tackle 
the question of correspondence between the almost exact 
phase space invariance present in these calculations and the 
persistence of density profiles in simulations of the signif- 
icantly more sophisticated cosmological context. We then 
outline the Fokker-Planck formulation used in this paper 
(Section and apply it to pure power law systems (Sec- 
tiontSjl. Section|S|discusses the evolution timescales on which 
the total mass distribution changes, relative to the dynami- 
cal friction timescale it takes for the clumps to sink in. These 
are then applied to the power law solutions in Section |7| 
In Section [H] the evolution of the individual components is 
considered in more detail. The consequences of the solu- 
tions obtained are outlined in several sections that follow. In 
particular, Section W7\ discusses the formation of a constant 
density core. This occurs in the background material if no 
stripping is invoked. The effect of stripping is dealt with in 
Section [T3l Both configurations correspond to approximate 
equilibria, with the ones involving stripping longer lived. 

Technical background material is outlined in several ap- 
pendices. Its inclusion was found necessary because of the 
variety of conventions, terminology and notation used in 
the literature of the Fokker-Planck equation as applied to 
gravitational systems. Moreover, to the author's knowledge, 
derivations of some of the results required for this paper have 
only been been published in the (French language) paper of 
Henon (1961). Alternative derivations are presented. 
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Figure 1. Density and velocity distributions of the total and background material composed of the lighter particles (denoted by DM). 
Passive evolution, refers to a system evolved from the same initial conditions, but with no massive particles. To reduce noise, the total 
distribution is averaged over three outputs, at 33.5, 32 and 30.5 dynamical times within the initial scale length r s . The background 
output, which does not contain high mass particles and is therefore less noisy, is given at 33.5 dynamical times, also at r s . Note that due 
to some initial evolution, resulting from the finite size of the system, this is slightly different from the scale length of the distribution 
shown in the figure. Since our comparison involves the passively evolved system, this effect is not of major importance. 



2 NUMERICAL EXPERIMENTS 

2.1 Invariance of the phase space distribution 

Fig. Q shows the density and velocity dispersion distribu- 
tions of systems of massive clumps and smooth background. 
The simulations are identical to one of the runs presented in 
El-Zant et al. (2004). The clumps are represented by soft- 
ened point particles whose 'size' is about 1/40 of the scale 
length of the Navarro, Frenk & White (1997) initial density 
distribution. There are 900 000 background ('light') parti- 
cles through which moves a system of 900 clumps with com- 
bined mass representing 20% the system's total. They are 
all started from the same initial density and velocity distri- 
butions. The system is spherically symmetric with isotropic 
velocities, and is evolved using a polar particle mesh code. 

As already noted by El-Zant et al (2004), the total den- 
sity distribution of the clump-background system remains 
virtually unchanged for many dynamical times (given in 
terms of the average density inside the initial halo scale 
length by td = 1/y G(p) s ); it is identical, except perhaps at 
the very centre and up to statistical noise, to the passively 
evolved system containing no clumps. This statement is also 
true at all times smaller than the timescale shown in the fig- 
ure (the invariance breaks down completely on timescales 
almost an order of magnitude larger, for reasons we discuss 
in the next subsection). 

We have tried different values of the NFW concentra- 
tion parameter c = r V i r /r s . The results presented in this 
figure refer to a system with c = 3.3. We ran, in addition, 
simulations with of c = 10 and c = 33.3, with almost iden- 
tical results. Different heavy particle masses and fractions 
were also investigated. As long as the problem is rescaled in 



terms of the relevant dynamical and dynamical-friction (see 
next subsection below) timescales, the situation described 
by Fig. remains for time intervals of the order of those 
shown. 

In El-Zant et al (2004) we had interpreted the clump 
distribution to be composed of baryonic material. In that 
case, the conclusions outlined in that paper still stand: the 
dark matter, made of lighter background particles, is heated 
and forms a constant density core; but the total distribution 
(interpreted to be composed of baryons and dark matter) re- 
mains constant. However, if the clumps are also made of dark 
matter, this invariance implies that the total distribution of 
dark matter remains constant. 

A rather remarkable result, transpiring from Fig. Q is 
that it is not only the physical density that remains un- 
changed, but also the velocity dispersion. The invariance of 
the velocity distribution is all the more remarkable because 
the initial 'temperature inversion' — the fact that the veloc- 
ity dispersion decreases towards the centre — is a thermo- 
dynamically unstable state of affairs, that should be washed 
out under the action of energy transfer; the system then 
evolving towards an isothermal distribution, as heat flows 
from hotter to colder regions. It so happens, nevertheless, 
that the 'heating' of the lighter particles is precisely balanced 
by a corresponding 'cooling' in the clump distribution, so as 
to keep the total 'temperature' constant. 

The constancy of the mass and velocity dispersion dis- 
tributions in a spherical isotropic system, where the distri- 
bution function is dependent only on energy, must be asso- 
ciated with an invariant phase space structure. In Section [3] 
we argue that this state of affairs, so stark in this idealised 
setting, may remain relevant, resulting in the observed per- 



© 0000 RAS, MNRAS 000, 000-000 



4 A. A. El- Z ant 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



r/r 

s 

Figure 2. Dynamical friction time for an NFW system with c = 
3.3, relative to the evolution timescale of the simulation of Fig. HI 
as a function of radius. It is assumed that the massive particles 
have 'size' l/40r s 



sistence of universal profiles, in the highly more complex cos- 
mological context. The remainder of the paper is an attempt 
to explain and derive the simulations presented here, and to 
quantitatively justify this basic contention. 



2.2 Parametrisation in terms of the dynamical 
friction time 

The state of affairs described by Fig.Qis metastable; in the 
sense it is long-lived, but not infinitely so. This should be 
obvious, since once the massive particles form a self grav- 
itating component, relaxation resulting from their mutual 
interactions will result in evolution. 

For timescales of order 200 rrj(r s ), we have verified that 
the center expands, as energy is transferred from the outer 
hotter regions. This is in accordance to what was found, for 
example, by Hayashi et al. (2003). Although we have not 
performed integrations for longer intervals, it is expected 
that it would then proceed towards core collapse. 

The notion of solid clumps cannot possibly be useful 
however in either of these regimes; the effect of stripping, 
which will transform clump material into background, en- 
sures that the former can never form an autonomous dy- 
namical system (cf., Section |3Jl 

We will see later in this paper (Section|S]and Section^, 
that the characteristic timescale for the evolution that sig- 
nificantly modifies the distribution of the (light and heavy) 
components, but keeps the total virtually invariant, is the 
local dynamical friction time. It will be useful therefore to 
parametrise the timescale of the simulation just analysed in 
terms of that quantity. 

To get a rough estimate of the relation between the time 
the system is evolved in Fig. Q and the average dynamical 



friction time within a given radius, we write the latter as 
(e.g., El-Zant, Shlosman & Hoffman 2001) 1 

E ri 

Ti5F(sim) = — ~ 0.75- — -t d , (1) 
E m A 

where r\ = M(< r)/m is of the total system mass inside 
radius r divided by the mass of a clump and A = ^ max is the 
Coulomb logarithm. We denote this timescale with the label 
'sim' to distinguish it from an analogous timescale td f that 
we derive later, in a more rigorous manner, from diffusion 
coefficients (Section |SJ. 

For a system with significant softening, as is the case 
here, 6 m i n is determined by the softening length. The max- 
imum impact parameter can be taken as the radius in- 
side which we are calculating the average dynamical friction 
timescale. We plot in Fig.|2]the dynamical friction time, rel- 
ative to the simulated timescale, for an NFW system with 
c = r v i T /r s — 3.33. 

Note that the dynamical friction timescale is smaller 
than the simulated one only in the very inner region, and 
then only by a factor of a few. (It increases again, because 
the softening dominates. Probably by then however, the for- 
mula above becomes entirely inadequate, the dynamics be- 
ing dominated by a non-Newtonian contributions to the soft- 
ened potential). Thus, in terms of dynamical friction time, 
the system is fairly young, despite it being of significant dy- 
namical age (parametrised in terms of td)- 

To explain the results of the simulation discussed in 
Section 12. II therefore, it suffices to show that the evolution 
timescales characterising the change in total phase space 
mass density of a system with central density cusp are much 
longer than the dynamical friction time. Noting also that 
most of the evolution in the lighter component in Fig. 
takes place deep inside the 1/r cusp, it should be possible, 
for the purpose of explaining the aforementioned results, to 
approximate the system as a scale free configuration. 



3 CONNECTING IDEALISED SIMULATIONS 
TO MORE COMPLEX SITUATIONS 

There are two glaring idealisations pertaining to the calcu- 
lations presented in the previous section: the heavier bodies 
are solid softened point particles; they also all have the same 
mass. The advantage, of course, is that the simplified models 
serve as to highlight and isolate the interesting phenomenon 
involving the invariance of the phase space distribution — 
and identify its origin in the dynamical friction interaction 
between a distinguishable background and a system of mas- 
sive objects, the number of which is large enough so that 
statistical reasoning may be relevant. But how do these toy 
models relate to the actual situation found in realistic simu- 
lations of halo evolution in the context of CDM cosmologies? 
In other words, would the discovered invariance survive, at 

Note that this equation approximates the density distribution 
as p ~ 1/r 2 . It also assumes that the velocity of the clumps 
is comparable to the dispersion of the background, which has a 
Maxwellian velocity distribution. 
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least in some approximate sense, the removal of the ideali- 
sations? 

In the following we attempt to make this connection 
The arguments are only suggestive; though most are broadly 
supported by the analysis based on the Fokker-Planck for- 
mulation presented in subsequent sections of the paper, fu- 
ture detailed numerical investigation should also be under- 
taken. 



3.1 Mass spectrum 

Subhaloes found in cosmological simulations do not all have 
the same mass. Indeed, their differential number distribution 
follows the rather steep mass function 4|r ~ m -9 / 5 (e.g., 
Gao et al 2004c). The mass is also a deciding factor as to 
how strongly a body will interact with a background — the 
rate of energy lost by that object varying as m 2 . It would 
thus appear that in the assumption of equal masses lies a 
drastic approximation, compromising the description of the 
evolution. 

Nevertheless, it may precisely be the steep mass spec- 
trum, and the sensitive dependence of the dynamical friction 
on the mass, that ensure that the situation described by the 
single (clump) mass model survives. For, by virtue of these 
properties, a separation of timescales will ensue; the heav- 
ier mass clumps sinking in first, while the lower mass ones 
are hardly affected by dynamical friction. If clump-clump 
interaction is not of major importance, this initial evolution 
will mimic that of the single mass model (clump-clump in- 
teraction will cause the dynamical evolution of small mass 
clumps to follow the background material). 

The total distribution of the system just described is 
also an equilibrium. Formally, we show in Section 2] that 
the fundamental equation deciding the presence of an equi- 
librium is linear with respect to the distribution of mas- 
sive objects: if a set of distributions, possibly corresponding 
to different mass species, is a zero of that equation, so is 
their sum (or any linear combination thereof). Thus, one can 
superpose an equilibrium distribution made of the evolved 
system, composed of the heaviest clumps plus background, 
and the distribution of the remaining (non-participating) 
clumps, which will remain unchanged if clump-clump in- 
teractions can be neglected (or, if such interactions are in- 
cluded, to various distributions of the remaining clumps). 

There is a caveat. Detailed invariance will require that 
the clumps representing substructure be initially distributed 
in a similar manner to the background. But the viability of 
a statistical description will be severely affected by the rela- 
tively small number of massive objects in each mass range. 
In particular, the most massive bodies, which initially con- 
tribute most, will be represented only by a few specimen. 
The invariance therefore would only hold in an ensemble av- 
eraged sense. Ma & Boylan-Kolchin (2004) found that sub- 
structure may, in individual runs increase or decrease the 
steepness of density profiles; and cosmological simulations 
do indeed show large scatter in the the profiles of identi- 
fied haloes (e.g., Klypin et al. 2001). More work, involving 
a series of simulations representing a statistical sample of 



initial sustructure distribution and concentrations, in the 
controlled context of isolated halo calculations, is in order. 

If the clumps come to completely dominate the mass 
distribution at small radii (as happens to the evolved sys- 
tem in Section r2.H . then the effect of dynamical friction will 
be limited to the most massive specimen — by the time, at 
any given energy, the dynamical friction coupling is strong 
enough to affect the next group of clumps, most of the back- 
ground particles that these clumps should couple with would 
have been removed. One is then left with a region that can 
be best described as a system of self gravitating clumps, 
with fluctuating and dissipative forces roughly balancing 
each other, but where energy transfer ('relaxation') leads to 
eventual evolution (cf. Sections IU1 and \TH . It is quite likely 
however that stripping would be effective in transforming 
much of the substructure into background, before such a 
situation transpires. This represents another form of cutoff 
for the dynamical friction-mediated coupling. 

Two question are clearly pertinent. They are concerned 
with the timescale for stripping and its effect on the total 
phase space mass distribution function. 

3.2 Stripping 

Satellites can no longer be regarded as solid clumps if any 
of the following processes are efficient 

(i) Particles are removed from the satellites due to prun- 
ing by the background (this includes the possibility of com- 
plete dissolution). 

(ii) Particles are removed due to weak encounters with 
other satellites. 

(iii) Satellites participate in strong mutual encounters or 
merge. 

The first of these mechanisms does not, to a very good 
approximation, change the total phase space mass distribu- 
tion function: particles are gently removed from a satellite; 
they leave with approximately the same velocity and posi- 
tion as the satellite. Thus, material with virtually the same 
location in phase space is transferred from clump to back- 
ground, the total phase space mass density remaining nearly 
constant. 

This will also be true of the second mechanism. Weak 
encounters will lead to the gentle removal of particles, they 
may also facilitate stripping by the main halo (cf. Penneru- 
bia & Benson 2004) , both processes through which the total 
phase space mass distribution is closely conserved. 

The situation is different however if violent encoun- 
ters and mergers are involved. It is possible to estimate 
the timescale associated with such events by averaging the 
clump distribution inside a given radius and invoking an ef- 
fective collision crossection. The envisioned 'collisions' will 
not necessarily lead to mergers (this will depend on their 
orbital parameters etc), but we will assume that when two 
satellites so interact, particles in their outer regions can 
be removed with significantly modified phase space coor- 
dinates.. 
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The collision mean free time for a randomly moving 
system of N spheres of radius d, confined within a spherical 
volume of radius r, can be written as (e.g., Saslaw 1985) 

1 r 2 

Tf =3Ndt Tc ' (2) 

where r c ~ 2td, is the time for an object to make one full 
crossing of the system. Clumps of different masses will have 
different diameters, but one can define an effective crossec- 
tion 



f d 2 dm 

J dm 

J; 



dm 



(3) 



Let the differential mass function ^ oc m 



9/5 , and the 

characteristic radius of a clump be d = cti B m 1 ' 3 , where the 
parameter B relates the virial radius and mass (it depends 
on the redshift and the cosmology; cf. Eq. Al of NFW) and 
Qi is the typical fraction of the virial radius of a clump that 
survives stripping by means of mechanisms i) and ii) in the 
list above. In this context we can write 

m~ 17 ' 15 dm 

(4) 



f 

D 2 2 Ji 

a — B «! 



f 



m -9/5dm 



2/3 



where we have assumed that the maximal mass in the clump 
distribution is at least a few orders of magnitude larger than 
the minimal one. 

1/3 

Confider a radius r — aBM v [ v , inside a host halo of 
virial mass Af v i r , and containing N(r) satellites. We can 
rewrite Eq @ in these terms, replacing d? by the effective 
crossection <r, 



Tf = 



N(r) 



«\ 2 / M vir \ 2/3 

cti) \m min / 



(5) 



If ai = a = 1, Tf w 20 T cr ~ 40 td (e.g., if N(r v i T ) is of order 
1000 and m m i n /M v i r w 10~ 6 ). This represents an average 
collisional timescale inside the virial radius. It is necessarily 
a lower limit, because it effectively assumes that stripping 
by the host halo, or via weak encounters, is negligible (hence 
ai = 1). 

As one averages over smaller radii (that is a < 1), our 
estimate of Tf will depend on how the ratio a/ai changes. 
A simple single power law that can be taken as a reasonable 
representation of the density distribution of host and satel- 
lites is p ~ 1 /r 2 (we discuss this further in Section|5Jl. In this 
case, it is easy to show (e.g., El-Zant, Kim & Kamionkowski 
2004) that the tidal radius of a satellite is proportional to 
its radial position inside the host halo — that is ai oc a. 
Since N(r) necessarily decreases with r, this suggests that 
the above estimate at the virial radius is indeed a lower limit 
— not only there, but at all radii. 

During a time interval r/, substructure continues to be 
radically modified due to stripping by the host halo; as can 
be seen from Fig. 11 of Taylor & Babul (2004) for exam- 
ple, where for almost all haloes, well over 90% of satellite 
mass is lost on such a timescale (note that the radial pe- 
riod used there P ra( j ~ 2t CI ~ 4td; and that the authors 
do not consider clump-clump interactions at all. Not even 
weak encounters). It therefore appears unlikely that strong 
encounters or mergers will be important at any stage during 



the evolution — haloes being affected by gentle stripping on 
a significantly shorter timescale. 

Finally, note that only a substructure satellite for which 
tdf < IOtd will be affected by dynamical friction before 
most of its mass is lost via stripping. From Eq. Q it should 
have ?7>0.01. Thus stripping eventually puts an effective 
end to the dynamical friction coupling. 

The discussion of the present section suggests that a 
situation whereas sinking clumps would come to dominate 
the mass distribution within a certain region (as in Fig. 
probably would not, in practice, arise — satellites having 
been significantly stripped long before such a configuration 
materialises — but that, nevertheless, the invariance found 
there will remain if 

(i) Stripping is due to gentle removal of particles from 
clumps due to their interaction with background or soft mu- 
tual encounters. 

(ii) Substructure distributions, properly modified to take 
into account the effect of stripping, also constitute (at least 
approximate) equilibria. 

The estimates of the mean free time presented here sug- 
gest that the first condition is satisfied. In Section 1131 we 
show that distributions satisfying condition 2 are indeed ex- 
cellent approximate equilibria. In fact, as may already be 
expected from the discussion above, such configurations are 
much longer lived, as they do not suffer from the problem of 
relaxation, characteristic of a system composed of self grav- 
itating clumps. 



4 FOKKER-PLANCK FORMULATION 
4.1 Background 

The Fokker-Planck equation has been discussed in the con- 
text of gravitational systems in several textbooks (Saslaw 
1985; Binney & Tremaine 1987; Sptizer 1987). The basic as- 
sumption is that the perturbations due to the granularity in 
a system are weak, local and random. The situation is then 
analogous to that of Brownian motion in a laboratory set- 
ting. The dynamical variables therefore undergo a random 
walk, the amplitude of which is limited by a deterministic 
force. 

Chandrasekhar (1943) reviewed these phenomena, both 
in the case of laboratory systems, where the approximations 
involved are much more easily justifiable, and the gravita- 
tional case; for which he derived an explicit formula for the 
deterministic force, that he dubbed dynamical friction. It 
proved remarkably successful in the face of numerical tests, 
despite the problems associated with justifying the assump- 
tions of locality and randomness (e.g. Zaritsky & White 
1988). 2 



2 Higher resolution simulations, testing the associated assump- 
tions in detail, are currently overdue nevertheless. For resonant 
effects, not included in the Fokker-Planck analysis, may mate- 
rialise at higher particle number (I thank Martin Weinberg for 
pointing this out). 
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The ratio of the amplitude of the fluctuating force caus- 
ing the random walk to that of dynamical friction is deter- 
mined by the mass of the object in question relative to that 
of the perturbers in the background. This implies that more 
massive particles will typically have smaller velocities; in the 
gravitational context, they will sink to the centre of the sys- 
tem. The timescale for this to happen we call the dynamical 
friction time. Eq. Q represents an estimate derived on the 
basis of the Chandrasekhar formula alluded to above. Other 
evaluations can be directly deduced from the Fokker-Planck 
formulation (Sec. |SJ. 

The energy deposited among the population of lighter 
particles causes their velocities to increase and their distri- 
bution to expand. The principal object of this paper is to 
show that these processes leave the total phase space mass 
density invariant for timescales significantly longer than the 
dynamical friction time (so as to be negligible on timescales 
of that order). In physical space, that would correspond to 
the expansion of the lighter particle system being almost 
exactly compensated by the inflow of massive clumps. 

Provided that the dynamical friction time is signifi- 
cantly longer than the orbital timescale (~ td), the Fokker- 
Planck equation can be averaged over particle trajectories, 
and expressed in terms of the integrals of motion in the 
smoothed out potential. The orbit averaged Fokker-Planck 
equation and its properties have been elucidated in detail by 
Henon (1961), and discussed in Binney & Tremaine (1987) 
and Spitzer (1987). Other works include that of Kohn (1979, 
1980) and Merritt (1983), the latter being particularly rele- 
vant to the situation at hand here. 

In the case when the configuration is spherical and has 
isotopic velocity distribution, as we will assume here, the 
six dimensional phase space is also spherically symmetric. 
The relevant integral of motion is the energy per unit mass 
E, which defines a radial coordinate in that space, in terms 
of which the distribution function and the Fokker-Planck 
equation may be written. 3 

We will not be using the Fokker Planck equation di- 
rectly here, but instead relations for the change in mass and 
energy, within a given surface E — const, that can be readily 
derived from it (e.g., Henon 1961 and appendices lAl and ITjl 
of the present paper). The advantage of this approach is 
twofolds. First, the Fokker-Planck equation itself allows for 
solutions with non-vanishing constant mass flux. These are 
steady state solutions, even though they are unphysical, un- 
less there is something to produce or absorb particles at 
the center of the system, a possibility we do not consider 
here. The second is that steady state solutions of the Fokker 
Planck equation do not guarantee thermodynamic stability, 
we will need to estimate the effect of energy transfer on the 
steady state (cf., Section f4. 41 . 



3 Note that E in this paper refers to the energy, and not to the 
binding energy, which has opposite sign, favoured in some of the 
aforementioned work. 



4.2 Mass change and the fundamental equation 
for the flux 

4-2.1 Mass change 

For the mass evolution we have (Appendix [XJ 



(6) 



dM(< E) dq 
d^ = ~ F + 9 m' 

where g is the phase space mass density distribution function 

— that is, the amount of mass contained in a volume element 
ArAv = (4ir) 2 r 2 v 2 dvdr is gArAv. This equation simply 
says that the rate of mass change inside the energy surface 
E is determined by the mass flux F — that is, the amount 
of mass, per unit time, that enters (negative sign) or exists 
the phase space volume bounded by the energy surface E 

— and the mass acquired or lost at the edge of this volume 
because of its variation with time. That phase space volume 
is related to the potential; since, in a spherical phase space, 



q{E) = (4jt)' 



2 2 , , 167T 

r v dvdr — — - — 



2 3 ■ 

r v dr, 



(7) 



and v — yj2{E — cf>) (here r max (£ ; ) is the maximum radius 
out to which a particle with specific energy E can travel). 
The 'density of states' p, the area in phase space of the 
surface with energy E, relates to the volume inside that 
surface by p = Accordingly, the mass enclosed within 
an energy interval [E, E + dE] is pgdE; and the total mass 
within the surface at E: 



(8) 



M(< E) = / pgdE 



Note that, even when the Fokker Plank flux F is zero, 
the mass distribution may be time dependent. This state 
of affairs is characteristic, for example, of smooth (collision- 
less) systems in configurations out of dynamical equilibrium. 
In this paper we assume that we are dealing not with the 
merger of a few clumps of comparable mass, but with a 
large number of them and a smooth background, with the 
combined system having reached a slowly evolving quasi- 
equilibrium state. Any evolution must therefore arise from 
the flux F describing the collisional evolution. We will there- 
fore ignore evolution driven solely by the second term on the 
right hand side of ©. 



4-2.2 General equation for the flux 

Comparison of the isotropic Fokker Planck equation in flux 
conservation form (Eq. IA6II and its standard format (e.g., 
Binney & Tremaine 1987; Spitzer 1987) yields 



from which follows the equation for the flux 

f=(d e -d ee ^) 9 , 

where 

D E =p(AE)-^p{(AEf) 



(9) 



(10) 



(11) 
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and 



D EE = - P ((AE) 2 ), 



(12) 



and where (AE) = Di and ((AE) 2 } = D2 are given by 
Eq. £3} and KJZl . 

The Fokker-Planck flux is due to the deterministic and 
random forces, acting on the components of the system due 
to their mutual interactions (as discussed in the last subsec- 
tion). The changes in energies that result from these forces 
are represented by the diffusion coefficients De and Dee 
respectively. It is possible to heuristically describe how their 
forms, and the accompanying expression for the flux, arise. 

When (11 LP and 1121 1 are inserted into llOL this expres- 
sion is seen to consist of three components. The first term, 
due to the component De g of HOL i.e. (AE) pg, refers to 
a shell in phase space of volume p (AE). On average, due 
to systematic change in energy (AE) , particles will traverse 
such a shell during a unit time interval. That is, particles 
crossing a surface at E = const in a unit time interval, b 
will come from this shell. The mass inside that shell is the 
volume multiplied by the phase space mass density. Hence 
the expression (AE) pg. 

The random fluctuations characterised by ((AE) 2 }, as 
opposed to their weighed gradients discussed below, can- 
not produce a particle flux unless the distribution function 
has a gradient — because particles will be as likely to enter 
or leave energy surface E. The RMS energy space 'length' 
traversed in unit time by a typical particle due to random 
fluctuation is S = ((AE) 2 } 1 / 2 . Consider a constant energy 
surface centered on the interval [E — S, E + 5]. The volume 
of each of the two shells bounded by E and surfaces at the 
extremities of this interval is p5; the difference between their 
average phase space mass density is <5^§; and so the differ- 
ence in mass they contain is pS x Since the walk in 
energy space described by ((AE) 2 ) is random, from each 
of those shells, on average, half the mass would transit in 
unit time through E (the other half crossing to still higher 
energies for the outer shell, or smaller energies in the case 
of the inner one). The magnitude of the resulting mass flux 
through E is then ±p((AE) 2 ) j& . It is directed outwards (in 
our convention positive) when the phase space mass den- 
sity is a decreasing function of energy — hence the term 
-D E e§§ in JT@. 

Systematic changes in energy can arise due to the de- 
terministic dynamical friction forces; or because the random 
fluctuations, weighed by the phase space available at energy 
E, p — p(E), are energy dependent. When this latter com- 
ponent is subtracted from total systematic energy change 
(AE)pg, the result represents flux whose sole source is dy- 
namical friction — that is, the term De g in llfty . 



4 Note that a flux will arise even if ((AE) 2 ) is energy indepen- 
dent; simply because, while particles will be as likely to cross en- 
ergy surface E from above or below, if p has a gradient, there will 
be more particles at energy intervals [E, E + dE] corresponding 
to larger p. 



4-2.3 Case of system with two components carrying 
masses m 2> fJ. 

In such circumstances, the average rates of change of the 
energy and its square are given by 1C3L 1C4L and 1C5II . It 
is then straightforward to show that De = for the light 
particles. For the heavy particles, it consists of the second 
term of Eq. (IU4t . 5 One can also see that Dee is the same 
for both species, and depends only on the distribution of the 
heavy masses. 

This is to be expected. Dynamical friction on the lighter 
particles is negligible; they can also be considered non- 
interacting, their evolution in energy space being determined 
by the random energy they, on average, gain from the mas- 
sive particles. For, if these heavier particles are removed, the 
remaining system, composed of light particles, is collision- 
less. 

The flux in light (background) particles is then 



Fu, = -D, 



dg^ 

■ dE ' 



(13) 



For the massive particles, on the other hand, we have 

F m = (De-Dee-^) g m . (14) 
Adding the two equations we get 

F = D E g m -DEE%. (15) 



By inserting the values of De and Dee into 1151 we 
obtain what we will call the fundamental equation for the 
flux 



—F 



pgdE+ \q g m dE - 
V J E 



,!E 1 f! (16) 



where we have assumed that the zero point of the potential 
is so chosen so that particles have positive energies (r = 
167r 2 G 2 lnA). 

The flux is a radial vector in an energy-parametrised 
spherical phase space. In the context of the conventions 
adopted here, positive flux, pointing outwards, corresponds 
to particle transport towards higher energies. 

Another equation that will be useful is the flux of back- 
ground material only. We have (by using [Tot 

which does not have a dynamical friction term, and thus 
invariably corresponds to mass outflow into higher energies 
levels. 



4.3 Linearity of the fundamental equation 

Now note that Eq. I i61) is linear in g m - That is, if a set of 
distribution functions g l mi (E) — mif l (E) are solutions, any 
function 



5 Note that our definition of De contains an extra mass factor 
as compared with that adopted by Merritt 1983, who prefers to 
multiply the term De in Eq. 11 01 by the mass. 
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g m = J2 A i9 l mi (E), (18) 
i 

is also a solution. Therefore one can construct solutions for 
the various components of a multimass system and add them 
up, assigning arbitrary weights, possibly derived from as- 
sumptions concerning the effects of mass segregation and 
stripping (cf. Section HS! . Thus our focus on two component 
systems does not lead to any fundamental limitation from 
the point of view of finding equilibrium solutions. 

This remarkable property is a derivative of the fact that 
the source term for the dynamical friction depends only on 
the total mass distribution; while the heating, though orig- 
inating only from the massive clumps, affects all particles. 
Note, nevertheless, that the evolution of components of dif- 
ferent mass will not be the same, since the dynamical friction 
flux for each individual heavy species i would be 

= ~A,m,g mi \ pgdE, (19) 
1 Jo 

while the flux due to the fluctuating force on members of 
that species is 

^ = Y^A imi [q J™ g l mi dE + j\g l mi dE^ ^.(20) 

Their ratio is Fr, F /F^ uc ~ m,/ mi. Heavier clumps will 
therefore lose energy to the lighter ones. 

From this it follows that the linearity of the fundamen- 
tal equation does not imply that the evolution of the system 
can be decomposed into the sum of the evolution of each 
species plus background separately! What it does imply is 
that, if each distribution corresponding to (e.g) distinct mass 
species happens to be an equilibrium solution for a given to- 
tal mass distribution, the combined system will also be in 
equilibrium. 



4.4 Energy change 

The amount of energy entering the energy surface E per unit 
time is 

dH{<E) — EE — f FdE + EF d 4 f F%E. (21) 
dt J dt J dt 

This equation is derived in Appendix^] (another derivation 
can be found in Henon 1961; cf. his Eq 2.27, 2.28 and 4.27) 6 . 

The first term of (I2H refers to the energy carried by 
the mass flux, while the second represents nonlocal contri- 
butions. The last two terms correspond to the variation in 
both these components due to the change in phase space 
volume, which results from change in the potential. They 
are analogous to the second term on the right hand side of 
Eq. @. 

6 Note that in Hcnon's paper S refers to our F and F refers 
to the distribution function. He also uses the term 'fundamental 
equation' to describe the Fokker-Planck equation instead of the 
flux equation, as the term is used here. 



This form takes into account the fact that particles en- 
tering into energy surface E from higher energy levels, at any 
given time, do not end up at the same energy Ej < E. Nei- 
ther do they all originate from the same initial energy level 
Ei > E. A distribution in Ei and Ef is always present (even 
if excessively large jumps are unlikely, because of the small 
deflection approximation assumed in deriving the diffusion 
coefficients) . The same goes for particles exiting energy sur- 
face E to higher energies. Thus, even in the absence of a 
net mass flux through a given surface E, energy can still be 
carried through it. This can lead to evolution. 

Because of the state of affairs just outlined, the evolu- 
tion driven by the energy flux will take the form of 'heat' 
transfer from 'hotter' regions of the system to those that 
are 'colder'. When a temperature gradient exists, statisti- 
cally, particles crossing E from above will end up at en- 
ergies which, though smaller than E, are larger than the 
compensating particles that crossed from below E to exit 
towards higher energies — again, even if the mass flux van- 
ished across E. This is why thermodynamic equilibrium, for 
a system of solid single mass objects, can only be isother- 
mal; and all open gravitational systems in virial equilibrium, 
so composed, are thermodynamically unstable (e.g., Saslaw, 
1985, 2000; Padmanabhan 1990; El-Zant 1998). Neverthe- 
less, we will show that in the case of the two component 
systems considered here, the evolution timescales associated 
with the energy flux can be quite long; and that, under cer- 
tain circumstances which may be associated with stripping, 
even exact equilibria that are not isothermal can exist. 



5 POWER LAW SOLUTIONS 

In this section we consider the case in which both the clump 
and the background distribution functions are power laws 
in the energy. The corresponding physical densities are also 
pure power laws. Perfect power law densities p ~ r 1 with 
— 1.5<7< — 1 are relevant representations in the inner re- 
gions of haloes, where the dynamical friction coupling is 
most efficient. Densities with indices —2 < 7 < — 1.5 are use- 
ful approximate forms at intermediate radii, up to the virial 
radius in haloes of small concentrations (as discussed in the 
penultimate and final paragraphs of Section 15.21 . We gen- 
erally have in mind a range — 2<7<0, although the case 
of 7 < —2, corresponding to the outer regions of haloes is 
briefly considered in Section [5.31 Strictly speaking the case 
of 7 = —2 requires separate treatment, because the potential 
is logarithmic in the radius. Nevertheless, the flux associated 
with profiles having 7^—2 continuously tends to the value 
at 7 = —2 (namely equilibrium). So will skip such separate 
treatment for the sake of brevity. 

5.1 General properties 

If the physical density can approximated by a power law 
such that 

p~r 7 , (22) 
the Poisson equation implies a potential is expressed as 
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Figure 3. Ratio of flux obtained from Eq. 1161 . to the first term 
in that equation, describing the flux of sinking massive clumps, 
for power law distributions following the forms given by Eq. 
and 1251 . The corresponding densities are p ~ r 7 for the total 
distribution and p m ~ r 7m for the massive clump distribution. 



■ = <pa r 



-0 



(23) 



with /3 = —(7 + 2). The distribution function that is an 
equilibrium of the unperturbed collisionless system takes the 
form (cf., Evans 1994; Evans & Collett 1997) 



go E 2 '^ 1 ' 2 . 



(24) 



Any distribution of a subset of particles having a phase 
space mass density function 

gm^gmo E 2/ ^- 1/2 (25) 
has corresponding physical space mass density given by 

p oc 

p m = iV2iT g m yjE - <f>(r)dE, (26) 



Pm{r) - 2/^ + 1/2- 
where the hypergeometric function 
Fi = F 1 {2/f3 m + 1/2, -1/2; 3/2 + 2//3 m ; 1) 



(27) 



(28) 



is a constant (that is with third argument equal to unity) 
when both g and g m are perfect power laws in the energy. 
The density of that component is then still a power law, 
with index 



7m 



P + 2 



(29) 



Finally, systems with scale free form for their total mass 
distribution function have density of states of the form 



dq 
dE 



■ p E 



1/2-3/ p 



(30) 



5.2 Mass flux 

A steady state implies that all time derivatives must vanish. 
In terms of Eq. JSJ this in turn requires that F — * 0. 
Substituting 12 It and 12511 into 1II6II one obtains 



F = Tmg g m o CE^-? + ?, 
where 



C = 



-p 



p- 



4-/3 



2/3„ 



PPr, 



3/3 - 6 I \4 + P 2/3/3 m - 3/3 m + 2/3 



(31) 



.(32) 



When P — /3 m Eq. I|31|l with F = reduces to the form 
studied by Evans & Collett (1997), thus their solution with 
/3 = —2/3 and 7 = —4/3 is relevant to clump-background 
systems considered here. As noted by these authors, the only 
other physical solution, that is a single power law, is the 
singular isothermal sphere (7 — > —2). 

When the assumption /3 = /3 m is relaxed, there is an- 
other solutions for each value of /3. It however corresponds to 
Pm > P- In the central regions (i.e, where E, r — > 0), where 
the dynamical friction coupling is strongest, and which must 
therefore be the focus of our analysis, the physical clump 
density would be larger than the total; which is unphysical. 

As an example, we plot in Fig.|H]the flux obtained from 
Eq. 131L normalised by the dynamical friction flux (i.e., by 
the first term of on the right hand side of that equation 
with C substituted in explicitly), for the case of /3 7^ f3 m 
with /3 = —2/3 — an exact solution at /3 = /3 m . In addition 
to this solution, there is another with /3 m = —1/2, that 
is 7m = —2. As E, r — > 0, this inevitably implies that 
the clump mass density is larger than the total (clump and 
background combined). This is clearly impossible (unless a 
low energy cutoff is introduced; see below). 

As one moves towards flatter clump density distribu- 
tions relative to the total, i.e. 7 m > —4/3, there is a net 
positive flux that increases sharply beyond "f m > — 1. Be- 
cause this flux is positive, and increases outwards, there is a 
decrease in phase space density; an expansion of the system. 
Moreover, this decrease in density is largest as one moves to 
smaller energy — suggesting that initial evolution proceeds 
in the direction of flatter density distribution. By inserting 
9m = 9 — 9n: in Eq. lib! , with g still given by 12 H and 
ff/j ~ g^oE 2 ^^" 1 ^ 2 , one can arrive at a complementary re- 
sult concerning the background component. This time, be- 
cause of the negative sign, one concludes that if the back- 
ground component has a power law distribution that is less 
steep than the total, there is compression of the total dis- 
tribution. These deductions are only suggestive; neither rig- 
orous nor comprehensive. They are nevertheless supported 
by the Monte Carlo experiments of El-Zant, Shlosman & 
Hoffman (2001), which show that systems with initially ho- 
mogeneous clump distributions expand. 

Predicting the initial direction of evolution in the case 
when 7 = 7,„ 7^ (—4/3, —2) is even less straightforward than 
the case when 7 = —4/3 and 7 m > —1 discussed above, 
because the flux diverges at small energies. It would appear 
that the particle distribution is rapidly depleted at these en- 
ergies, destroying the cusp and forming a core; while the dis- 
tribution steepens outside that (expanding) region (because 
dF/dE < 0). In practice, a small core in the spatial clump 
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distribution must be always assumed to exist, because sub- 
haloes have finite size. It corresponds to a low energy cutoff 
in the phase space distribution function (Section which 
would ensure well behaved forms for the flux at E — > 0. Such 
a small energy cutoff in the clump distribution does not lead 
to appreciable deterioration in the accuracy of the solutions 
— as we show in Section Hol m connection with the issue of 
stripping. 

An important feature of Fig. [3] is that although for 
f3 = /3 m the only exact solution (other than the singular 
isothermal sphere) is that with exponent 7 = —4/3, for 
7 < — 1 the normalised flux is quite modest, suggesting that 
these configurations may be long lived. In order to make such 
statements more precise, we consider, in Section^ the ques- 
tion of evolution timescales. They determine how small the 
flux needs to be in order for a given approximate solution 
to be considered relevant for the timescales of interest. 

Power law densities with 7< — 3/2 are steeper than 
inner cusps seen in numerical simulations. They can nev- 
ertheless be retained as good represntations to the density 
profiles of cosmological haloes over a large radial range. For 
example, as can be seen from Fig. 1 of El-Zant & Shlosman 
(2002), fits to cosmological halo density profiles correspond 
to softened 1/r 2 distributions over more than three orders 
of magnitudes in the softening length. These profiles do not 
incorporate faithful reproductions of the phase space struc- 
ture of the inner cusp, since they lack the required temper- 
ature inversion in velocity space, but they should represent 
excellent approximations at intermediate scales (on which 
velocity profiles are approximately isothermal), up to the 
virial radius for haloes of small concentrations. 

Thus, the existence of approximate solutions in the 
range —2 < 7< — 1 may be interpreted to mean that 
gravitational systems with density profiles similar to those 
of cosmological haloes are approximately invariant under 
substructure- induced interaction — since dynamical friction 
coupling at r ^> r s should be negligible, except for the most 
massive satellites. Nevertheless, as we note below, if r S> r s , 
and power laws with 7 > — 2 no longer represent reasonable 
approximations to the density distribution, the normalised 
flux can be large, even divergent. This, in principle, may 
have an effect if dynamical friction coupling on some satel- 
lites is not entirely negligible, and if a significant fraction of 
the host halo mass is contained in these regions. 




-2.4 -2.3 -2.2 -2.1 

Y 

Figure 4. Same as in Fig. El but for 7 < -2 (Note that Ea.ll6l 
has to be slightly modified in this case, as explained in text). The 
flux diverges for 7 < —2.5. 

of the blow up is the third integral in 1161 . When the total 
mass distribution is a power law, so that q ~ p E, it can 
be seen to correspond (in absolute value) to a mass weighed 
average of the energy of the clump distribution (recall that 
AI(E) = pg). The excess of low energy clumps causes their 
average energy to diverge for 7 < —2.5. The result is a 
divergent positive flux. For 7 < —3 the first integral also 
diverges (because the mass does). 

In practice, a cutoff is introduced, because 7 increases 
as r — > r s . The severity of the net flux will then depend on 
the mass inside the region where 7 falls significantly below 
— 2. Cosmological haloes do not generally contain a large 
fraction of their mass in regions beyond a few scale length 
r s . The dynamical friction coupling is also exceedingly small 
in these low density regions. Nevertheless, the possibility 
that evolution on a fraction of a dynamical friction time 
occurs in the outer regions of highly concentrated haloes 
cannot be completely ruled out within the context of the 
present analysis; neither can the direction of any prospective 
evolution be determined. 



5.3 Distributions with 7 < 2 

Although power laws with indices 7 ~ — 2 may be useful 
representations of the density up to the virial radii for cos- 
mological haloes of small concentration, beyond a few scale 
length r s , 7 — > —3. 

For 7 < —2, it is appropriate to replace the lower 
bounds in the first and third integral in Eq. I16H by — 00, and 
the upper bound in the second integral by 0. Also E — > — E 
in equations J^J, and (3DJ|. 

The normalised flux is shown in Fig. 0] in the range 
—2.5 < 7 < —2. Only in the immediate neighbourhood of 
the singular isothermal sphere is the flux comparable to the 
case of — 2<7< — 1. It diverges for 7 < —2.5. The source 



5.4 Energy flux 

A thermodynamically stable steady state requires that the 
time derivatives in Eq. 12H vanish. This in turn implies that 
(EF — J FdE) — > 0. For power law systems this requirement 
translates to 

D = C f 1 " + a ) =° ( 33 ) 

(where C is given by I32H . For f3 — (3 m , this is impossi- 
ble, except for j3 —* —2, corresponding to (an unphysical) 
completely homogeneous system, or f3 — + (the isothermal 
sphere). All other solutions are thermodynamically unsta- 
ble. To determine whether these are nevertheless long lived 
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for the timescales of interest we must derive associated evo- 
lution times. This is done in the next section. 

But, for 13 7^ p m , there are exact solutions that are 
thermodynamically stable for all values of j3, provided that 



P m = 



6/3 
4-/T 



(34) 



That is, the massive particle distribution is slightly less steep 
than that of the total. This cannot correspond to any evolu- 
tion with solid clumps — since, in this case, the clumps 
sink in, replacing the background to dominate the mass 
distribution near the centre of the system (r, E — > 0). It 
may be relevant however in the presence of stripping (Sec- 
tion EUi a process whereby clump material is transformed 
into background. Note that, because these configurations 
have p m m (3, the corresponding mass flux is also quite small. 



6 EVOLUTION TIMESCALES 

6.1 Timescales associated with mass flux 

A standard procedure in physics is to derive a characteristic 
timescale by freezing the physical parameters of a system 
at a given moment and estimating the evolution time from 
there. In our case, using Eq. (JSJ, a timescale for the evolution 
of the total mass distribution because of the existence of a 
non-vanishing mass flux can then be denned as 



TF = 



M(< E) 
\F{E)\ ■ 



(35) 



which is the time interval the flux at energy E takes to 
significantly modify the total mass distribution inside that 
energy surface. 

Suppose that, within that energy surface E, the mass 
distribution is initially dominated by the smooth back- 
ground component. Then the fluctuations in the potential 
resulting in random heating (described by the coefficient 
Dee) have negligible effect on the clumps. They therefore 
sink in under the action of dynamical friction on a timescale 
similar to that given by Q. A closely related quantity, the 
time on which the mass in sinking clumps significantly mod- 
ify the original mass in clumps inside energy level E is given 
by 



M m (< E) 
\F DF \ ' 



(36) 



where Fdf as usual corresponds to the first term or the right 
hand side of Eq. 1161 . 

The dynamical friction time itself can be arrived at by 
using Eq. HI IP . again neglecting self interaction terms arising 
from clump-clump heating due to their mutual encounters 
(i.e., the term involving {(AE) 2 )). One gets 



E 



TDF 



(AE) 



pE 
~De~ 



pEg„ 
\Fdf 



But M m (E) = p(E)g m (E); so we have 

EM m (E) 
TDF = ^fV 



(37) 



(38) 



Since M(< E) is generally smaller than EM(E), this 
timescale is larger than that obtained from 1361 : it involves 
the effect of the flux on the local mass distribution of clumps, 
instead of its influence on the distribution inside the phase 
space surface bounded by E. 

One can also define an analogous timescale for the evo- 
lution of the total mass distribution. The ratio of these 
timescales would be 

M m (E) \F\ 



R = tdf/t - 



(39) 



M(E) \F DF \' 

It determines how rapidly the total mass evolution takes 
place relative to the dynamical friction timescale; which de- 
scribes how fast the distribution of one of the components, 
namely the clumps, changes. 

Suppose now that the initial system consists of single 
power law configurations (i.e, in terms of the terminology 
of Section 7 m = 7). The dynamical friction flux for this 
initial system then is 

-mV „ 



Fdf = 



c, 1/3+1/2 



1-1//3 

(obtained by insert ing |24l and l25l into [T^t : in terms of which 
one gets an explicit form of the dynamical friction timescale: 



T DF (t = 0) = 



1/p 



E 



1/2-2//3 



(41) 



mTgo 

We omit the explicit reference to the initial time in what 
follows, even thought this is implied. 

For such a system to keep its total phase space mass 
density g nearly constant (Eq. l24l remaining a good approx- 
imation) , for timescales of the order of the initial dynamical 
friction time, it is necessary that the absolute value of the 
dimensionless flux 

tdf = g m (E) F g m0 (E) F 
r ~ g(E) F DF go(E) \F DF \' 1 ' 

be significantly smaller than unity. Note that this same re- 
sult can be inferred from 13511 and 13611 . 

Note also that, although this equation assumes a power 
law form for the total distribution, the flux F can correspond 
to arbitrary distribution for any of the individual (i.e, back- 
ground and clump) components, which are necessarily time 
dependent. 

The quantity F e defines a relative error, quantifying 
the validity of approximate solutions. More precisely, a so- 
lution is a valid approximation, that is the system's total 
phase space distribution stays approximately invariant for a 
timescale t if \{F e ) t \^^ < 1 for all relevant energies. At 
any given moment one can define a dimensionless character- 
istic evolution timescale 



Te = l/\F e \; 



(43) 



it expresses the number of dynamical friction times the to- 
tal phase space mass distribution can remain approximately 
invariant. A system can be considered unstable under the ac- 
tion of the coupling induced by the presence of substructure 
clumps if r e w 1 — since, in that case, the total mass distri- 
bution changes if the configuration is evolved long enough 
for dynamical friction, the principal manifestation of such 
interaction, to be effective. 
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As an example, the system studied in S ect ion \'2 . 1 1 must 
have (F e ) >C 1 for energies corresponding to radii that are 
evolved over a dynamical friction time or longer. 

Fig. Q also suggests that the timescale for each of the 
components to be significantly modified is comparable to 
the dynamical friction time. This is consistent with the re- 
sults in this section. Note however that the replacement of 
background by clumps does not proceed linearly in time — 
in which case it would take a time ~ -jj 8 - tdf{E) for the 
clumps to dominate inside some energy surface E. The some- 
what more detailed analysis of Section[S]shows that the evo- 
lution is actually exponential in tdf (at least initially). This 
explains why the population exchange is virtually complete 
on timescales of the order of a dynamical friction time, even 
if the mass fraction in clumps is quite small. 

Finally, note that if the mass fraction in heavy clumps 
is not negligible, then the effective dynamical friction flux 
decreases by a factor of g^o/go, because part of that flux 
is compensated by heating due to clump-clump interaction 
(c.f., Eq. 1160 . The corresponding dynamical friction time 
increases by an inverse of this factor. The effective evolution 
flux can then be defined as 



Fc.rtf — 



ffmO 



F 



ff M o \Fdf\ 



(44) 



Unless otherwise stated, we will have in mind a system where 
the background dominates, so that F e<e fi « F e , and will 
be using Eq. 1421 ; since it corresponds to a definition of 
the dynamical friction time that is closer to conventional 
usage based on the Chandrasekhar formula, which considers 
only one massive particle in an infinite background medium 
(Eq. 0is based on that formulation); obviously, conversion 
between F e and F e B s is trivial. 



Eq. 1431 is necessarily large compared to the dynamical fric- 
tion time. 

For example, for the mass ratio of the simulations of 
Sec. 12.11 g m o/go = 0.2, and so F e ~ 2% in the cen- 
tral region, where the evolution takes place and where 
7 ~ — 1. This means that when both mass components 
are distributed in power law density configurations with 
—2 < 7 < — 1, the total phase space mass distribution would 
evolve over ~ 50 dynamical friction timescales, very large 
compared to timescales of interest (at r s , this corresponds 
to w 15000 td(j's); which is far larger than the age of cos- 
mological haloes). 

When f3 = /3 m the energy evolution flux can be written 

as 



F e E = 



g m a 



P 2 (2 + P) 



g (1-/3X1 -20)08 + 4)' 



(47) 



This reaches a maximum of 0.111 for 7 = — 1. It is 

so ' 

smaller for other values of > 7 > —2. Scale free systems in 
this range of 7 are therefore also long lived from the energy 
transfer viewpoint — that is, the time for energy transfer 
to affect the total distribution is much longer than the dy- 
namical friction time, (note that, for 7 = — 1, F e — F e E, 
although there appears to be no reason to suppose this not 
to be entirely coincidental). 

The mass evolution flux is not negligible when 7 > — 1 
— that is, for flat initial density distribution — especially 
if the mass fraction in clumps is not small. In this case one 
must use Eq. (1441 . with g m o ~ g^o; which, for 7> — 1, 
implies that -F e , e s ~ 1- Therefore, a system which starts 
with a weak central cusp, and a significant fraction of its 
mass in clumps, is unstable if the coupling due to dynamical 
friction is non- negligible. The same conclusion applies to the 
evolution driven by the energy flux when 7 ~ — 1. 



6.2 Timescale associated with energy flux 

The same attitude adopted in the subsection above allows 
one to define timescales associated with the energy flux. For 
example, the expression analogous to 1351 is 



TE 



EM(< E) 



EF- f FdE\ ' 
One can also derive a corresponding 'evolution flux' 
gmo EF — j FdE 



FeE = 



.9i) 



E\F DF \ 



(45) 



(46) 



7 EVOLUTION TIMESCALES OF 

APPROXIMATE SCALE FREE SOLUTIONS 

We now apply the results just derived to power law density 
distributions with "f m = 7 (considered in Section^. 

The values of F/\Fdf\ in the case of these scale free so- 
lutions are shown in Fig- 01 F° r ^ 2 < 7 < — 1 the absolute 
value of this ratio is < 0.111. For models with single power 
laws g m o /go < 1 determines the mass density of clumps rel- 
ative to the total mass density. Eq. 1421 then implies that 
evolution timescale associated with the mass flux (given by 



8 SEGREGATION INSTABILITY 

Initial scale free solutions can only be adopted as zeroth 
order approximations. As the system evolves, the coupling 
between the two components will cause the light background 
particles to move out of lower energy levels and be replaced 
by the massive clumps. In addition, the latter can lose mass 
to the background under the action of stripping. The goal of 
the remainder of this paper is to show that the total phase 
space mass density can nevertheless remain nearly constant, 
approximately keeping the power law form, even when the 
individual distributions are modified by the aforementioned 
processes. 

Eq. @ , for the background particles alone, can be writ- 
ten as 



c>M M (< E) 
dt 



-F.{E)+g^ t , 



(48) 



where _F M is given by Eq. 1171 . Suppose now that the total 
phase space distribution is largely unchanged over timescales 
of interest; that is, F e <C 1, and so g continues to be ex- 
pressed in terms of its initial value form, given by Eq. 1241 . 
The second term on the right hand side of 1481 is then also 
negligible. 
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We assume a generalised power law solution for the 
background # M = g^o(t, E)E 2/l3 ~ 1/2 . This, in itself, is not 
an approximation, since the function g^o is arbitrary. How- 
ever, in order to obtain a solution, we assume, in addition, 
that for the purpose of evaluating the integrals in (1171 , g m 
keeps its original scale free form, given by 1251 . For the first 
integral in Eq. 1171 this is an excellent approximation at any 
given time t and energy E such that t/rDF(E) < 1; since, in 
that case, the effect of the dynamical friction coupling is just 
becoming important at energy E, and is therefore unimpor- 
tant at larger energies (recall that tdf ~ 7J 1/ ' 2 ~ 2/ ' ,a ). 

For the second integral, one notes that when the total 
mass distribution remains in scale free form q — pE/ (3/2 — 
3/(3). This integral then represents a mass weighed average 
of the energy in massive clumps (recall that M(E) = pg). 
It will decrease as the the latter lose energy because of dy- 
namical friction. Nevertheless, if this occurs principally at 
low energies, where tdf is smaller, the change is not dras- 
tic. More formally, integration by parts transforms it into 



7/(3/2 - 3/(3) = E / pg m dE 



dE 



pg m dE, (49) 



where the scale free form for (Eq. 1301 remains relevant, in 
accordance with our assumption that the total distribution 
is practically time independent. Using Eq. JSJ, we get 



7/(3/2 - 3/(3) =M m (<E)E- / M m (<E)dE. (50) 

Jo 

Mass conservation implies that the first term in this 
expression is constant if the effect of dynamical friction is 
small inside E. It is also always larger than the second 



M m (<E) E > / M m (< E) dE, 



(51) 



where the equality corresponds to the case when all the mass 
is concentrated near E — — rather unlikely given the as- 
sumption that, at energy E, the effect of dynamical friction 
is just being felt at some energy E ^ 0, (it would imply a 
sharp discontinuity in the M m (E) profile). Moreover, 7 is 
always smaller than the first integral in Eq. 1171 . which as 
we have seen is virtually constant. In the initial scale free 
system its contribution to the total outward flux is a frac- 
tion (2//3 + l/2)/(2//3 + 1/2 - 2+ 1/(3), which is less than 
1/3 for (3 > — 1 (i.e., for 7 < —1). This ratio can only be- 
come smaller on timescales for which the inflow of massive 
particles across E does not significantly change their total 
mass inside 75, but where redistribution inside this surface 
renders their phase space density profile steeper (again this 
follows from Eg. 1501 , 

The use of the initial scale free form for g m in evaluating 
the integrals on the right hand side of (1171 is thus justified 
when t/r df(E) < 1. Since the total mass inside energy sur- 
face E is still unaffected on that timescale, A7 M (< E) can 
also be be computed using the initial scale free distribution. 
In this context, Eq. 1481 transforms into 

^=-A(g - 9lJ , ) giM0 , (52) 
where 



A = ml a(!3) E 2/ ?~ 1/2 , 
with 

a((3) = (1 - 1/(3) 



4-/3 

2/0 + 1/2 2-1/(3) \ 3(3-6 



(53) 



(54) 



a parameter of order unity. Eq. (I.">2li can be integrated to 
give 



= X (E)e~ soAt . 

<?m0 



(55) 



Here we have assumed, as usual, g m = g — g^ (and so g m o — 
go — .9mo)- The function \ is m principle arbitrary, but is in 
fact independent of energy when the initial ratio of light to 
heavy particles has this characteristic. And so we have 



<?m0 



ffmO 

The timescale 
1 



-90-4* 



(56) 



(57) 



determines the interval on which the exchange of popula- 
tions, associated with the massive clump inflow accompa- 
nied by expulsion of background, takes place. This segrega- 
tion timescale is of the order of the dynamical friction time. 
More explicitly, using Eq. (EJl, 

I^=a((3)(l-l/(3). (58) 



9 EFFECT OF LOW ENERGY CUTOFF ON 
THE ACCURACY OF EQUILIBRIUM 
SOLUTIONS 

9.1 Behavior at low and high energies 

The results of the previous section suggest that dynamical 
friction coupling causes depletion in the distribution of back- 
ground particles; first at the smallest energies, and then at 
progressively higher energy. 

Suppose that there is some cutoff energy E c below 
which no background particles exist. At energies E <C E c , 
and still assuming that the clumps are indestructible (no 
merging or striping), there are again equilibrium solutions 
which are single power laws — this time for a one compo- 
nent system composed of the clumps, since one can now 
rewrite the fundamental equation only in terms of the 
clumps (by putting g — g m in Eq. 1161 . The dimension- 
less evolution fluxes (Eq. 1421 and 1461 are now significantly 
larger, since g m o/go — > 1; the associated timescales corre- 
spondingly smaller. Nevertheless, this state will still evolve 
on a timescale that is > 10 tdf, provided 7< — 1. Only 
when 7 is significantly larger, corresponding to flat density 
profiles, will the evolution timescale be comparable to the 
dynamical friction time. 

It is important to note here that the ultimate effect of 
dynamical friction by a smooth background on solid clumps 
is to produce a self gravitating core made of these clumps, 
with the total density distribution largely unchanged — and 
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not the sinking of the clumps to the centre as may be ex- 
pected had the background been of infinite mass. No matter 
how small the fraction of mass in clumps initially, the ratio 
will increase exponentially on a dynamical friction timescale 
(as we have shown in the previous section), with the clumps 
eventually coming to dominate inside some (time depen- 
dent) energy level E — E c . This puts an effective end to 
effects induced by dynamical friction from the background, 
which has been expelled from that region, and to the sinking 
of clumps to the physical centre of the system. 

For E E c the system would still retain its original 
distributions of clump and background particles, the dynam- 
ical friction time being too long for there to be any effect. 
The question therefore is what transpires at the boundary; 
the transition region between the two equilibria. 



9.2 Boundary region: sharp versus gradual cutoff 

9.2.1 Mass flux 

Suppose that g M = for E < E c , but remains unmodified 
for larger energies. 7 By writing g m = g — g M in Eq. II (it 
one can separate it into two components: one involving only 
functions of g and another involving both g and g^. If it is 
assumed that the total phase space mass distribution func- 
tion g remains a perfect power law, then the flux due to 
the one species term (involving only g), can be made arbi- 
trarily small — by choosing an appropriate value for (i.e, 
~ —2/3 or — > 0). The two species component (involving 
g and g^) has no such solutions. In fact, for E < E c there is 
the residual flux 



Ri_ 

mV 



q I g^dE 



dg 
dE 



P b(J3) E\ 



2/g+i/z E -i/0 



where P = Po go ff M o and 
4-/3 



(59) 



(60) 



(3/3-6)(2//3 + l/2)' 

which is always nonzero. 

Dividing by the dynamical friction flux, we obtain the 
associated evolution flux (cf., Eg. 1421 



F e {Ri) = 



9n0 



H0) (E c \ 
1 - 1/3 \ E ) 



2//3+1/2 



(61) 



go 1-1/8 

As may be expected from the discussion of the preceding 
subsection, for E <C E c this residual flux becomes small; 
the system of clumps, thus represented, is near equilibrium. 
For E ~ E c (and ~ —g^o/g ~ — 1) however, it is smaller 
than unity but not negligible. 

For E > E c there is another non-zero residual term 
resulting from the abrupt cutoff. It is given by 

7 Of course, strictly speaking, these conditions cannot be re- 
alised simultaneously, since the total mass in light particles should 
be conserved. However the phase space volume increases quite 
steeply with energy (~ E 3 / 2- so that the migration of low- 
E particles has little effect on the distribution there (as in Fig. HI 
where drastic changes in the inner distributions of light particles 
has little effect on the outer regions where they moved out). In 
what follows we will be assuming this approximation holds. 



Ro _ 
mT ~ 

with 
c(3) = 



< r , .11 )!<L=Poc(f3)Et 1/fl E^-^(62) 



4-/3 



(3j9 - 6)(2 - ■ 
The associated evolution flux this time is 
g,o c(j3)_fE c \*-^ 

IP 



F e (Ro) = (f ) 

go 1-1/0 V E I 



(63) 



(64) 



It has a sign opposite to that of F(Ri) and decreases away 
from E — E c significantly faster. It is also smaller at E ~ E c , 
but still is not completely negligible. 

Both residual terms can be made far smaller if one re- 
places the abrupt cutoff with a gradual one. We do this by 
introducing a variable cutoff energy E c = E C (E). This con- 
ceptually amounts to dividing the background material into 
a continuum of populations, each with its own cutoff energy, 
such that g^o{E c ) changes in an interval dE c by an amount 
dg^o- We will suppose that at some energy E c = E™ ln , 
<?mO = 0; and at another, E c = _E™ ax , it takes its initial value. 
The evolution fluxes can then be written (by summing over 
all the populations) as 



F e (Ri) = 
and 

F e (Ro) = 



b(0) I 



9f,0(Bmax) £,2/0 + 1/2 



9 M 0(B) 



dg^o 



1 - 1/0 ff0 E 2 /5+l/2 



c(0) jr {E) E!-^dg, 



1-1/0 goE 2 - 1 /! 3 



(65) 



(66) 



We can also rewrite Eq. (16511 and 16611 in terms of the 
dimensionless variable X = E c /E. For any phase space point 
with energy E the 'inner flux' from points with E c > E is 
then given by 



F e (Ri) = 



b(0) 1 
1 - 1/0 go J 1 



X 2 /0+ i/2 rdj^ dXj (67) 



where X max = E?**/E and X min = Ef™/E. There will 
also be an 'outer flux' from points with E c < E. It is 



F e (R 



o 



1 - 1/0 g 



-1//3 / dgu.0 \ 
\ OX J 



dX. 



(68) 



Note that, both these integrals are minimised when g^o is 
flat, that is (-§^-j «C 1, when g^o is largest. By choos- 
ing functions with decreasing logarithmic derivatives, one 
can make the integrals as small as one likes. The fact that 
the two terms are always nonzero also means that they tend 
as to cancel each other. As we will see below, it is easy to 
find functional forms that render the total evolution fluxes 
quite small. 



9.2.2 Energy flux 

In this case the relevant quantity is EF — J FdE, which 
replaces F in calculations analogous to those just described. 
It then follows ('using EfH and I59|l that 
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Figure 5. Residual evolution fluxes (calculated using |rj71 and l68t due to exponential cutoff in the phase space mass distribution function 
described by Eq 1711 . Left panel: /3 = —2/3 (corresponding to physical density p ~ 1/r 4 / 3 ). Right panel: /3 = —1 (p ~ 1/r). When 
j3 = — 1 we take into account that even systems without the low energy cutoff are only in approximate equilibrium (hence the term 
Fe{g ~ ff/j))- The total phase space mass distribution function g is assumed to keep its original scale free form (Eg. 1241 . 



F Ee (Ri 



1-/3 



F e {Ri 



and, in an analogous manner (this time using Irl2l 



F eE (Rc 



2/0 - 1/2 



F e (R c 



(69) 



(70) 



Thus, for /3 ~ — 1, the inner residual energy flux is smaller, 
and the outer larger, compared to the corresponding mass 
fluxes by only a modest factor. Since we will be interested in 
showing that the fluxes are actually smaller by an order of 
magnitude, or more, than those that would lead to signifi- 
cant change of the total phase space mass distributions, over 
timescales of interest, these factors are of secondary impor- 
tance. In the illustrations that follow we will be focussing 
on the mass flux. 



10 ACCURACY OF EXPONENTIAL 
SOLUTIONS FOR THE INITIAL 
EVOLUTION 

Equation l|56|l was arrived at under the assumption that 
the total phase space mass density remains constant during 
that evolution; that is r e 3> r s . We now check this. Self 
consistency requires that when we insert back the solution 
into till , the evolution flux is indeed small. 

We suppose that each element of the system can be 
evolved according to Eq. 156H . According to the discussion 
of Section |H1 this is a good approximation at points where 
the effect of dynamical friction is just becoming significant. 
Nevertheless, the distribution at smaller energies, where the 
dynamical friction time is smaller, need not follow this form. 
This affects the calculation of the third integral in 1161 and 
therefore involves a further approximation, even at points 
where t<TDF (note that it does not affect first integral 



which depends only on the total distribution). But since, 
as also discussed in the aforementioned section, this integral 
(denoted by I) is insensitive to the precise distribution in- 
side energy surfaces E for which t/rDF(E) < 1, and in any 
case does not dominate the flux, this approximation should 
be valid when calculating that flux for points where the dy- 
namical friction is just becoming important. 

Eq. 1561 can be solved for g m o at any given energy E 
and time t to give 



5m0 = 



+ (ML) e-t/r. 



-t/T 



9o- 



(71) 



Substituting this form for g^o, using Eq. I57H and 15311 . 
m lISTIl and |HH1 and ( since this functional form does not 
admit sharp cutoffs) letting X m i n — > and X max — > oo, one 
obtains the residual fluxes. These are shown in Fig. |S] in 
terms of 



i 1 ) 



(72) 



and for (g m0 /g ) t=0 = 0.2. 

Since the exponential evolution was deduced under the 
assumption that material beyond our reference point at 
X = 1 was largely unaffected, they are expected to be highly 
accurate only for s<2. For this is the value of s that cor- 
responds to one dynamical friction time for j3 ~ — 1 (from 
Eg. 1551. 

The results of Fig. do indeed show that for s < 2, 
the average value of the evolution timescale (the inverse 
of the evolution flux shown) correspond to ~ 20 tbf - 
which means that, up to one dynamical friction time, the 
timescale for the total mass to evolve away from its ini- 
tial phase space distribution is twenty times longer. The ap- 
proximation breaks down completely when t/roF \{F e )t \ ~ 
l/2s |(-F e )s| ~ 1, where the angular brackets refer to aver- 
ages of F e over times smaller than t. From Fig. |S] it may 
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Figure 6. Evolution fluxes (calculated via Eq. 1421 , describing the accuracy of approximate solutions of the form 1731 , for (3 = —2/3 
(corresponding to p ~ 1/r 4 ' 3 ,) shown on the left panel, and (3 = —1 (p ~ 1/r) (shown on the right). The total phase space mass 
distribution function g is assumed to keep its original scale free form (Eg. 1241 , 



be estimated that the exponential solutions are relevant for 
the whole period of evolution that is plotted. They are a 
good approximation (~ l/2s |(^ ? e)s| ^$0.2) up to s ~ 4, cor- 
responding to about two dynamical friction times. 

Note that, when s is very large, at X ~ 1 the system 
tends to the one-component equilibrium state described in 
Section r9.ll Therefore, the accuracy of the exponential solu- 
tion does not continue to deteriorate as s increases. However, 
because the gradient of g^o is quite large, the decay in the 
inner flux approximately follows the sharp cutoff form given 
by Eq. 1611 . and the tendency towards the one component 
solution is slow. But the form given by 1711 has been derived 
under the assumption that s ~ 1, and is therefore no longer 
relevant. In the following, we present approximate solutions 
that describe a smooth transition between the two equilibria 
discussed in Section r9.ll with corresponding evolution fluxes 
that are very small at all points. 



11 SMOOTH SOLID CLUMP SOLUTIONS 
VALID FOR ~ 10 - 100 t df 

Since, for power law forms of the total mass distribution 
function g, there only corresponds two solutions for the in- 
dividual components that can be also expressed in terms of 
powers of the energy (Section Q>J, no solutions can be writ- 
ten as an infinite power series in the energy, and be exact. 
Nevertheless, as we show here, some such solutions can be 
exceptionally excellent approximations. 

For example, we have tried solutions of the form 8 



8 We have also tried other algebraic (e.g., ~ (1 + (E / Eq)™) -1 ) 
and smooth exponential forms (e.g., ~ (1 — e -kE/E c ^ £ or g y 
with similar results. 



^=^)(l- (1 J w ), (73) 

where Eq is some scaling parameter and g IJ ,o{E 00 ) corre- 
sponds to the unperturbed value (i.e. to <? M o(i = 0)). They 
represent systems where the background material is de- 
pleted at small energies due to interaction with solid massive 
clumps. 

We insert g m = g — into Eq. (1161 . with g, the to- 
tal phase space mass density assumed to remain constant. 
Consistency requires that the evolution flux (Eq. I42H is 
small. This substitution breaks the equation into two ad- 
ditive terms. In the case when /3 — —2/3, the g — g terms 
have exactly vanishing total flux. For /3 = — 1 there is an 
evolution flux of 0.111. It will add to any evolution flux re- 
sulting from the g — g M term. 

The total evolution fluxes are shown in Fig. |S| (again for 
the case where (g m0 /g ) t = — 0.2). As may be expected from 
the discussion in Section [9. II when (3 = —2/3, they tend to 
zero as E/Eo <C 1 (corresponding to the one component 
equilibrium state); and as E/Eo 3> 1 (where the initial, two 
component equilibrium, still holds) . As may also be expected 
(this time from Section [§J, the absolute values of the fluxes 
peak at the transition between these two regimes, and are 
larger when the transition is most rapid. Note nevertheless 
that Eq evolves on the local dynamical friction timescale, 
as the instability front moves out to higher energies on the 
segregation time scale t 3 of Section|S| Therefore no given re- 
gion in energy space would remain, during the whole period 
of evolution, in the regime where the peaks lie. Like in the 
previous section where we discussed the initial exponential 
evolution, it is the average value of the flux that is relevant. 

When (3 = — 1, and E/Eo <C 1, the residual flux con- 
nected to the g — g (one component) term causes the total 
flux to tend to 0.111. The associated evolution timescale is 
of the order of ten dynamical friction times. The energy flux 
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also becomes important on such timescales. Significant evo- 
lution in the total mass distribution is therefore expected. 
We have verified by means of simulations (like those of Sec- 
tion l2.f \ that evolution in the total distribution does indeed 
occur on these timescales. ft is initially in the direction of 
expanding core (as found by Hayashi et al. 2003). We have 
not integrated our systems for longer times, but core col- 
lapse, driven by the evolution of the (by then completely) 
self gravitating system of clumps, probably occurs. Both 
these regimes however are well beyond any situation where 
the representation of substructure in terms of solid clumps 
can be expected to be of any use. For stripping would have 
radically modified the distribution and internal structure of 
satellite haloes. 



12 THE FORMATION OF A CONSTANT 
DENSITY CORE 

The existence of a low energy cutoff in the phase space dis- 
tribution will entail the removal of some particles of the 
affected species from the central region of the system in 
physical space; it does not however imply the removal of all 
particles. Instead, the radial mass density distribution of the 
background will develop a nearly constant density core. For, 
as can be seen by replacing the lower integration limit by E c 
in Eq. (1271 . there will always be a radius where E c 3> <j>(r). 
For power law potentials of the form <j> — (feor -13 , the transi- 
tion radius is r t = (f>o/E c . Deep inside this region the density 

I- a „E !/H1,! 

tends to the constant value p c = 4V27T 2/ ,3+1/2 — ■ 

This is of course what is found for the background dis- 
tribution in Fig. (left panel), where the density of back- 
ground particles does indeed flatten off at small radii. In El- 
Zant, Shlosman & Hoffman (2001) and El-Zant et al. (2004) 
we had interpreted the clumpy component to be made of 
baryons; and assumed that stripping was negligible and that 
this state of affairs actually materialised in practice (the ra- 



tionale being that baryonic material has formed by dissipa- 
tion and is therefore dense enough to survive). 

In contrast, if the clumps are made of dark matter, there 
will be continuous transformation between the two compo- 
nents, as stripping transforms the clump material into back- 
ground. Cosmological simulations actually suggest that this 
process may be so effective that the cutoff in the energy 
distribution in practice probably occurs in the clump distri- 
bution, and not the smooth background. For the density dis- 
tribution of the clumps consistently shows evolution towards 
a flattened core (e.g., Ghigna et al. 2000, Fig. 10). 



13 STRIPPING 

As discussed in Section f3.2l stripping does not significantly 
change the phase space mass density — because stripped 
particles are likely to escape from the clumps with very 
small kinetic energy and near zero binding energy (relative 
to the clump), the amount of mass in the phase space energy 
range [E, E+dE] remaining nearly constant (the mass in the 
clumps in that energy range decreases, but is compensated 
by an equal increase in background). 

Therefore, if the interaction due to dynamical friction 
leaves the system in equilibrium, it will remain so in the 
presence of tidal stripping. One nevertheless still needs to 
show that satellite distributions that are affected by strip- 
ping also correspond to long lived equilibrium states for the 
total mass distribution. 

Suppose now the stripping turns most of the mass in 
clumps into background materials for energies E<Eq. Then 
the role of the two components are reversed from what was 
described in Section HT1 Now, as appears to be the situation 
in evolved cosmological systems, it is the satellite distribu- 
tion which would have a lower energy cutoff and develop a 
core (in physical space) , while the background would have a 
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cusp. In that case, instead of substituting g m = g — p M into 
Eq. 116L as we have done above, we can directly put 



9mo(E) 



(7mo(-Eoo) 

(1 + E/E a y 



(74) 



Since this simply represents a direct exchange of the role 
of the two components compared to what is described by 
Eq. , we expect the situation to be similar — except that 
the evolution flux is now expected to be even smaller, be- 
cause the principal source for evolution, the massive clump 
distribution, vanishes at small energies, where the coupling 
is strongest. As can be seen from Fig. |7J this is indeed the 
case. 

Since there is now virtually no clump-clump interaction 
at small energy, there is no longer a relatively large flux as 
E — *• when = — 1. The energy transfer flux is also always 
small. Physically this simply means that, if stripping is ef- 
ficient over timescales smaller than the segregation time t s 
(Eq. 1571 . the evolved regions tends towards a non-evolving 
collisionless state as Eo increases, instead of the self gravi- 
tating system of satellites expected to transpire if the strip- 
ping time is arbitrarily long (i.e., when the clumps can be 
considered solid). 

Of course, in a realistic system, the distribution of 
clumps would be more complicated. There will be a mass 
spectrum, and each species may have a different distribu- 
tion. Nonetheless, because of the linearity of the fundamen- 
tal equation 11611 with respect to g m (cf. Section 14.31 . one 
can add any number of approximate solutions of the form 
1741 (or, for that matter. 1751 and still end up with a system 
that is close equilibrium — the errors, that is the evolution 
fluxes, will simply add linearly. 



14 CONCLUSION 
14.1 Principal results 

The principal conclusion of this paper is that the phase space 
mass distribution associated with density distributions akin 
to the 'universal' halo profiles are approximately invariant 
under the action of the interaction induced by the presence 
of substructure satellites; which is shown to lead to negligible 
evolution over timescales comparable to halo lifetimes. This 
conclusion applies to the central region ( p ~ 1/V), where 
the dynamical friction coupling is strongest, and up to radii 
where the profile can be approximated as p ~ 1/r 2 . 

Since haloes in cosmological simulations build up from 
smaller components; are continually incorporating infalling 
material; and retain significant substructure throughout 
their evolution, this is a necessary condition for the observed 
universal forms, apparently present at all redshifts and mass 
ranges — even if the initial collapse naturally leads to the 
'right' distribution, evolution driven by the presence of sub- 
structure will modify it, unless the total mass distribution 
is shown to be invariant over relevant timescales. 

The physical mechanism behind our claim is simple. 
Clumps and background are coupled via dynamical friction 
interaction. The clumps move in, the lighter background par- 
ticles move out, the total remains the same. This equation, 



which we have shown in terms of distribution functions in 
energy space, naturally leads to invariance in physical space. 
It is unmodified by the effect of inter-substructure interac- 
tion in the form of weak encounters. 

We presented simple simulations which showed that for 
dozens of dynamical times (inside the initial scale length 
r s ), the physical as well as velocity space total distributions 
remain invariant under the dynamical friction mediated in- 
teraction between a (solid) clumpy component and a back- 
ground of far lighter particles — despite the fact that, during 
that same timescale, the distribution of each component is 
radically modified (section l2.H . 

The rest of the paper had the two principal goals: 1) to 
explain the results of these simple simulations and 2) To dis- 
cuss, both from the qualitative and quantitative viewpoints, 
the relevance of such idealised models to the more complex 
context of cosmological simulations. 

The conceptual and computational framework em- 
ployed is based on a Fokker-Planck formulation, dividing 
self gravitating structures into two components: a back- 
ground made of light particles and a system of clumps rep- 
resenting substructure satellites. The latter are assumed to 
interact with the background via dynamical friction, and 
among themselves through weak gravitational encounters. 
Although we usually have in mind a system where the back- 
ground dominates, this is not an a priori assumption of the 
model. 

It shows (Section |SJ that there are exact scale free so- 
lutions for cusps with both components having scale free 
density distributions p ~ r 7 , and 7 = —4/3. This solution, 
already found for single mass systems by Evans & Collett 
(1997), is generalised for the case of clump-background sys- 
tems considered here. As long as the massive clumps are 
much heavier than the background particles, any mass dis- 
tribution of the former is allowed; by virtue of the linearity 
of the fundamental equation for the mass flux (Sec 14.31 . 

Power law distributions, in the whole range of 
— 1 > 7 > — 2, are in approximate equilibrium, both in terms 
of the mass and energy transfer; the total mass distribution 
in phase space remaining constant for many dynamical fric- 
tion times (Sec. 0. Solutions with 7 ~ —2 do not correspond 
to any power law found in the central regions of cosmologi- 
cal haloes. But they are faithful representations of the phase 
space structure at intermediate radii (up to the virial radius 
for haloes of small concentration; Section Nevertheless, 
when r>r s and 7 > — 2 power laws can no longer be taken 
as faithful represntations of the density distribution, there 
can, in principle, be a large mass flux towards higher en- 
ergies, if dynamical friction is not completely negligible at 
these radii, and if they contain a significant fraction of the 
halo mass (Section 15.31 . 

Even when the total mass distribution function is con- 
stant, there is continuous evolution in the distribution of 
each of the components. Two cases may be distinguished, de- 
pending on how efficient the stripping of the massive clumps 
is. 

In the case when the clumps are considered as solid ob- 
jects, we show that the evolution of each of the components 
away from the initial power law function takes the form of 
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an exponential instability, which we termed the 'segrega- 
tion instability', and which takes place on a characteristic 
timescale comparable to the local dynamical friction time. 
It involves the replacement of lighter particles, which gain 
energy and exit lower energy levels, by the massive clumps 
(Section 0. 

This results in a situation whereby the distribution of 
light background particles develops a low energy cutoff in 
energy space — with a corresponding constant density core 
in physical space (Section 1121 . in line with the evolution of 
the idealised iV-body models discussed in Section \'2. II Long 
lived solutions for such models are given in Section [Tl"l Situ- 
ations whereas substructure is formed by dissipationless col- 
lapse are unlikely to be represented by such models for any 
significant period of time. Nevertheless, these may be rele- 
vant if the clumps are made of significantly denser baryonic 
material; which is more resistant to dissolution via strip- 
ping, because it dissipated during collapse (e.g. Gao et al. 
2004b). One then recovers the situation described by El-Zant 
et al. (2004), where a core developed in the light dark matter 
particles but the total density distribution remained largely 
unchanged. Similar conclusions were reached by Gao et al 
(2004a) in a much more sophisticated cosmological setting. 

If the timescale for mass loss from clumps — stripping 
— is significantly smaller than the segregation time, a cut- 
off in the satellite distribution function may result instead, 
with an accompanying constant density core in the physical 
space distribution of that component. This appears to be 
what happens in the context of dissipationless cosmological 
simulations. 

In Section EH we had argued that stripping does not 
modify the total mass distribution function. Therefore, if 
the dynamical friction interaction between clumps and back- 
ground, and weak encounters between clumps, do not modify 
that distribution function, it suffices to show that solutions 
that take into account stripping — that is ones incorporat- 
ing a low energy cutoff in the clump distribution — can be 
long lived. Such solutions, valid up to a thousand dynami- 
cal friction times (much longer than the age of cosmological 
haloes), are given in Section [T3l 

Finally, it is interesting to note that, whereas in the 
case of a single mass system, or a multimass systems with- 
out stripping, the only solutions with no energy transfer 
(i.e., thermodynamically stable configuration), are isother- 
mal spheres, the case may be different when there is a mass 
spectrum and stripping is effective (Section 15.41 . 

14.2 Concluding comments 

Several issues raised in this paper would seem to require 
numerical work, either by means of the Fokker-Planck equa- 
tion or direct simulations, to be settled in a satisfactory 
manner. The two approaches are complementary: the full 
simulations including fewer physical assumptions and the 
Fokker-Planck models corresponding to more controlled cal- 
culations that are computationally far less intensive. The 
question concerning statistical effects arising from the fi- 
nite number of subhaloes, and the contention that stripping 
does not modify the total mass distribution function (both 



discussed in Section [3J, need to be empirically verified by 
means of full simulations, involving live TV-body satellites. 
Other issues are best treated by a combination of the two 
approaches. 

The effects of a mass spectrum and stripping may be 
explicitly included in Fokker-Planck models (Merritt 1983). 
It is also possible to incorporate a prescription for the evolu- 
tion of the halo mass distribution — as was done for exam- 
ple by Nusser & Sheth (1999). These authors used a stable 
clustering formulation, supplemented by an algorithm for 
the effect of dynamical friction on substracture and associ- 
ated back reaction on the main halo. Because they have as- 
sumed that energy lost by sinking satellites is redistributed 
homogeneously among parent halo particles, their conclu- 
sion was that sinking substructure causes a steepening of the 
cusp, invalidating the 'universality' assumption. Their basic 
model nevertheless remains relevant and can be coupled to 
the Fokker-Planck formulation, which provides a much bet- 
ter representation of the effect of energy deposition by sub- 
structure (which in fact is better approximated as a local 
phenomenon) . 

Of prime interest is the issue of stability, touched upon 
in the introduction — whether the effectively dissipational 
interaction due to the presence of substructure may lead to 
an 'attractor' like behaviour, with systems tending to pre- 
ferred configurations similar to those observed in the large 
simulations. There are several related points that are raised 
by the material presented here. 

The fundamental conclusion of this paper, the imparl- 
ance of the total phase space mass density distribution under 
the action of the dynamical friction coupling, strictly speak- 
ing holds only when both components initially have the same 
phase space distributions. Simple considerations involving 
the direction of phase space fluxes and their gradients (Sec- 
tion^ suggest that power law systems may expand if the the 
clumps distribution rises less steeply than the background, 
and contract if the situation is reversed. What is available in 
terms of numerical data tends to confirm this. For example, 
El-Zant, Hoffman & Shlosman (2001) ran Monte Carlo sim- 
ulations where solid substructure was given homogeneous 
spatial initial conditions inside an NFW halo; the combined 
system was started from virial equilibrium and left to re- 
lax towards detailed dynamical equilibrium; the dynamical 
friction, modelled on the Chandrasekhar formula, was then 
turned on. The total density distribution decreased. Heuris- 
tically, one may explain this from the fact that the binding 
energy available, per unit mass, is larger the more diffuse dis- 
tribution. Nevertheless, this issue is of sufficient importance 
so as to merit further work, involving a series of simulations 
whereby the spatial distribution of substructure within the 
halo, as well as the internal structure of individual clumps 
is varied. 

The material in Sections |S] and |7| suggest that sin- 
gle power law systems with flat cusp (7 > — 1) are only 
marginally unstable — they would evolve on a dynamical 
friction time only when 7-^0 (which may be taken to 
represent the nearly homogeneous centre of a system) and 
when the mass fraction in clumps is comparable to that in 
the background. Whether any evolution actually occurs will 
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clearly depend on if stripping is effective enough in order to 
quench the evolution, by modifying the mass ratio and cut- 
ting the source of the flux before any actual evolution takes 
place. It is of interest to determine under what (if any) con- 
dition evolution does occur, and whether this happens in the 
direction of the cusps characterisitc of cosmological haloes. 

Finally, there is the issue of whether the large evolution 
flux associated with the very outer parts of haloes (7 < 
—2; see Section f5.3|l has any significance in determining the 
structure of the outer profiles. Dealing with these regions 
may also require the incorporation of velocity anisotropies. 
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APPENDIX A: EQUATION FOR MASS 
CHANGE 

In general, a kinetic equation with a collisional term can be 
written as 



dg _ (dg\ 

dt vat A 



(Al) 



where g is mass density distribution function — the amount 
of mass inside a unit volume element in a six dimensional 
phase space. When the collisional term on the right hand 
side is zero, one recovers the collisionless Bolzmann equa- 
tion: the dynamics conseerves the phase space mass density 
along the motion. 

The collisional term can be written as 

'dg- 



Vat A 



(A2) 



where F p is the mass flux through phase space, defined as 
the amount of mass crossing unit area per unit time, and 
V p is the phase space gradient. This is simply a statment of 
conservation of mass in phase space — the change in mass 
density inside a phase space volume due to the effect of 
encounters is the a difference in the amount of mass entering 
and leaving as a result of these encounters. 

For a system where the distribution function depends 
only on energy, that is one whose phase space distribution 
is spherically symmetric, one can rewrite l|A2(l as 



v at v con 



1 OF 
pdE' 



(A3) 



Here F is the mass flux travesing the energy surface E. Since 
this does not have unit area, one divides by the phase space 
area of this surface, hence the factor 1/p. 

The term on the left hand side of (IAH can be written 
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1-g-V.V* (A4, 

which is the collisionless Boltzman equation written in terms 
of the phase space velocity V of a particle. When the phase 

ldq. 

p at ' 



space is spherical that velocity is E = — ; so that we have 



dg_ 

dt 



dg_ _ 1 dq dg 
dt pdtdE' 



(A5) 



Inserting this into 1A1L taking into account IA3I) . one gets 

(A6) 



dg 
dt 



1 OF 1 dg dg 
~pdE + p~dtdE' 



This can be rewritten, in terms of M(E) — pg, as 



dM(E) 
dt 



1 dF d 
pdE + dE 



which, upon integration, gives 
dM{< E) 



dt 



-(F(E)-F(0)) + F(t) 



d_ 

dE 



(•S) 



(A7) 



(A8) 



For an isolated system in quasi-equilibrium the arbi- 
trary flux F(t) = 0. If there are no source or sinks of mass 
at the centre of the system then also F(0) = 0. This is the as- 
sumption we adopt here. We nevertheless note that, strictly 
speaking, this cannot not be true of the approximate scale 
free solutions presented in Section since, for these cases, 
the flux diverges as E — > 0. One therefore has to assume a 
small core that breaks the scale free solutions at the centre. 
In practice this is naturally realised when the finite size of 
the clumps is taken into account. Under the zero central flux 
assumption we have 



dM{< E) 
dt 



-F(E) + g 



dq 
dt' 



(A9) 



APPENDIX B: 
CHANGE 



EQUATION FOR ENERGY 



The rate of energy change of a population of particles in 
an interval [E, E + dE] can be related to the rate of mass 
change inside that interval by 



dH{E) _, dM(E) 
— hi 



(Bl) 



dt dt ' 

or, in terms of the total mass and energy inside surface E, 



d_ dH{< E) 
dE dt 



E 



d dM{< E) 



dE 



dt 



(B2) 



Integration (by parts on the right hand side and setting ar- 
bitrary time dependent functions and central fluxes to zero) 



gives 



dH{< E) 
dt 



E 



dM{< E) 
dt 



dM(< E) 
dt 



dE, 



(B3) 



which, when used in conjunction with <A9t yields the re- 
quired energy change. 

Note that, in moving from the local energy H(E) to the 
integrated one H{< E), we have neglected the interaction 
potential of the particles, that is we have simply added their 



energy. This is consistent with the context of the Fokker- 
Planck formulation, since the assumption of locality in re- 
lation to the encounters producing the energy changes (c.f. 
Section already incorporates the idea that the evolution 
is cuased by encounters which proceed independenly of the 
form of the mean field, which is responsible for this inter- 
action. This is no longer true if singificant energy changes 
are caused by an evolving mean field. However, for our pur- 
poses, this point is not of central importance; the rationale 
being that we assume that the system is in quasiquilibrium 
and evolves only due to the encounters. The source of any 
evolution is then the Fokker-Planck flux. 



APPENDIX C: DIFFUSION COEFFICIENTS 

In general, diffusion coefficients represent the averge rates 
that a certain quantity changes over time. In case of the 
isotropic orbit averaged Fokker-Planck equation they rep- 
resent the average rate of change of powers of the specific 
energy E, averaged all orbits with energies in the interval 
[E, E + dE]. (Note that, as always with this equation, this 
involves the assumption that the actual rate is such that the 
change during an orbital time is small, thus permiting the 
averaging of quantities over orbits). 

In the weak encounter (Fokker-Planck) approximation, 
two diffusion coeficients are relavant: D\ = (AE) and 
D2 = {(AE) 2 ). For a test particle of mass m t affected by 
encounters with field particles of mass mj and distribution 
function //, they are given by (Spitzer 1987; Theuns 1996) 



£>i = Fmj 
and 



f f dE-^L q f f dE 
e P m f Jo 



Lh - 2 \:,„i ( i j q /. di: ■ Z J f f dE 



(Cl) 



(C2) 



where T = 16-7T G In A, A being the Coulomb logarithm; 
q is given by J7J and p = Jfe (note that in this paper we 
are assuming that the zero of the potential is taken so that 
energies are always positive). 

In a system consisting of several mass species, field par- 
ticles will include all other particles of the same species, 
plus those of the other species present. The form of the 
diffusion coeficients is additive in the distribution of these 
different species, it is thus easy to generalise relations ob- 
tained for a two component system. We therefore consider 
such a system; with constituents having masses m and fi. 
The phase space mass distribution functions are g m = mfm 
and g M = We will assume that the phase space mass 
densities of the two components are comparable; that is g m 
and g M are roughly of the same order, while at the same time 
m 2> fi. For the average energy change of the light particles 
one therefore has 



(AE„) = -Ym I g m dE, 

J E 

while for the massive species it is 



(G3) 
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(AE m )=Tm(^J g m dE-jJ pgdE^j, (C4) 

where g = g m + g»- 

The average square changes in specific energy are 

((AS) 2 ) = 2 Tm (- [ qg m dE+l [ g m oIe) . (C5) 

\P JO P J E J 

This relation holds for both the light and massive species. 
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